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Abstract Inspired by the theoretical results on optimal preconditioning stated by Ng, 
R.Chan, and Tang in the framework of Reflective boundary conditions (BCs), in this 
paper we present analogous results for Anti-Reflective BCs, where an additional tech- 
nical difficulty is represented by the non orthogonal character of the Anti-Reflective 
transform and indeed the technique of Ng, R.Chan, and Tang can not be used. Nev- 
ertheless, in both cases, the optimal preconditioner is the blurring matrix associated 
to the symmetrized Point Spread Function (PSF). The geometrical idea on which 
our proof is based is very simple and general, so it may be useful in the future to 
prove theoretical results for new proposed boundary conditions. Computational re- 
sults show that the preconditioning strategy is effective and it is able to give rise to a 
meaningful acceleration both for slightly and highly non-symmetric PSFs. 
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1 Introduction 

Image deblurring problems H fTTlfT2l[T3l represents an important and deeply studied 
example belonging to the wide area of inverse problems. In its simplest form, the 
deblurring problem consists in finding the true image of an unknown object, having 
only the detected image, which is affected by blur and noise. In this paper we deal 
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with the classical image restoration problem of blurred and noisy images in the case 
of a space invariant blurring: under such assumption, the image formation process is 
modelled according to the following integral equation with space invariant kernel 

g(x) = J h(x-x)f(x)dx + r](x), xeR 2 , (1.1) 

where / denotes the true physical object to be restored, g is the recorded blurred and 
noisy image, 77 takes into account unknown errors in the collected data, e.g. mea- 
surement errors and noise. We consider a standard 2D generalization of the rectangle 
quadrature formula on an equispaced grid, ordered row-wise from the top-left corner 
to the bottom-right one, for the discretization of (JTTTJi. Thus, we obtain the relations 



in which an infinite and a shift-invariant matrix = [hi-j](i,j)=((i l ,i 2 ),{ji J2))' *' e ' : 
two-level Toeplitz matrix, is involved. 



Though ( 1 .2 1 relies in an infinite summation since the true image scene does not 
have a finite boundary, the data g, are clearly collected only at a finite number of 
values, so representing only a finite region of such an infinite scene. The blurring 
operator typically shows also a finite support, so that it is completely described by a 
Point Spread Function (PSF) mask such as 



IPSF 



\hh ,h\u =-a, a, .(i=-oi ai (1-3) 



where i 2 > for any i\ , 22 and Yj=-q h — 1, i= (h,h), <l— (<7i;<?2) and the normal- 



ization is according to a suitable conservation law. Therefore, relations ( 1.2 1 imply 

1 

gi = E h sfi-s + f]u h = l,.-.,n\,h = l,...,m, (1.4) 

s=—q 

where the range of collected data identifies the so called Field of View (FOV). 

As in \\.2\ , we are assuming that all the involved data in ( 1.5 1 are reshaped in a 
row-wise ordering, so that the arising linear system is 

Af = g-t] (1.5) 

where A £ ^N(n)xN(n+2g) j s a fj n j te principal sub-matrix of A ro , with main diagonal 
containing hop, f G R N ( n+2q \ g, r\ € M. N ( n > and with N(m) = m\ni2, for any two-index 
m = {m\.m.'i). 

According to ( |1.4} , the problem is undetermined since the number of unknowns 
involved in the convolution exceeds the number of recorded data. Thus, Boundary 
conditions (BCs) are introduced to artificially describe the scene outside the FOV: the 
values of unknowns outside the FOV are fixed or are defined as linear combinations 
of the unknowns inside the FOV. In such a way ( |1.5[ ) is reduced to a square linear 
system 

A n f = g-ri (1.6) 
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withA„ G R N ^ xN ("\n = (n 1 ,n 2 ),N(n)=mn 2 and/,g,7] G R Ni "\ 
Though the choice of the BCs does not affect the global spectral behavior of the 
matrix, it may have a valuable impact both with respect to the accuracy of the restored 
image (especially close to the boundaries where ringing effects can appear) and to 
the computational costs for recovering / from the blurred datum, with or without 
noise. Notice also that, typically, the matrix A is very ill-conditioned and there is a 
significant intersection between the subspace related to small eigen/singular values 
and the high frequency subspace. 

The paper is organized as follows. In Section [2] we underline the importance of 
boundary conditions and we summarize the structural and spectral properties of ma- 
trices arising in the case of Reflective and Anti-Reflective BCs. Section|3]is devoted 
to the presentation of theoretical results relative to the explicit construction of the 
optimal preconditioner for the restoration problem with Anti-Reflective BCs. In Sec- 
tion [4] we report computational results with respect to two deblurring problems, the 
former having a slightly non-symmetric PSF and the latter having an highly non- 
symmetric PSF, using the proposed preconditioning (combined with Tikhonov filter- 
ing) for Landweber method. Finally, some conclusions and perspectives are drown in 
Section |5] 

2 The role of boundary conditions in the restoration problem 

Hereafter we summarize the relevant properties of two recently proposed type of BCs, 
i.e., the Reflective ifTTl and Anti-Reflective BCs |20|, with respect both to structural 
and spectral properties of the resulting matrices. Indeed, the use of classical periodic 
BCs enforces a circulant structure and the spectral decomposition can be computed 
efficiently with the fast Fourier transform (FFT) |6|, but these computational facil- 
ities are coupled with significant ringing effects [4|, whenever a significant discon- 
tinuity is introduced into the image. Thus, the target is to obtain the best possible 
approximation properties, keeping unaltered the fact that the arising matrix shows an 
exploitable structure. Reflective and Anti-Reflective BCs more carefully describe the 
scene outside the FOV and give rise to exploitable structures. Clearly, several other 
methods deal with this topic in the literature, e.g. local mean value [19| or extrapo- 
lation techniques (see 1 14 1 and references therein). Nevertheless, the penalty of their 
good approximation properties could lie in a linear algebra problem more difficult 
to cope with. Hereafter, as more natural in the applications, we will use a two-index 
notation: we denote by F = \fi, ,J , , the true image inside the FOV 

and by G = \gh,h]h=\,..., ni ,h=\,...,n 2 the recorded image. 

2. 1 Reflective boundary conditions 

In El Ng et al. analyze the use of Reflective BCs, both from the model and the 
linear algebra point of view. The improvement with respect to Periodic BCs amounts 
to the preserved continuity of the image. In reality, the scene outside the FOV is 
assumed to be a reflection of the scene inside the FOV. For example, with a boundary 
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at x\ — and x 2 = the reflective condition is given by f(±xi,±X2) = f(x\,X2). 
More precisely, along the borders, the BCs impose 

fh,l-h = fh,h> fhm+h = fhm+l-hf * an y '1 = '2 = l,---,<?2 

and, at the corners, for every ij = 1,... z'2 = 1, . . . ,#2 the use of BCs leads to 

/l— 12 = fh,i%i fn l +i l M2+i2 = fni + l— i'i,n2+l— '2' 

f\—i\,n%+ii fi\.n2+\— (2 > /ni+ii,l— 12 i'i,!2' 

i.e., a double reflection, first with respect to one axis and after with respect to the 
other, no matter about the order. 

As a consequence the rectangular matrix A is reduced to a square Toeplitz-plus- 
Hankel block matrix with Toeplitz-plus-Hankel blocks, i.e., A„ shows the two-level 
Toeplitz-plus-Hankel structure. Moreover, if the blurring operator satisfies the strong 
symmetry condition, i.e., it is symmetric with respect to each direction, formally 

h^i=hi for any i = — q, . . . ,q. (2.1) 

then the matrix A„ belongs to DCT-III matrix algebra and its spectral decomposition 
can be computed very efficiently using the fast discrete cosine transform (DCT-III) 
EH . More in detail, let %, = {A„ G R N( -" )xN ("\n = (m,n 2 ),N(n) = n x n 2 \ A„ = 
R„A n Rl} be the two-level DCT-III matrix algebra, i.e., the algebra of matrices that 
are simultaneously diagonalized by the orthogonal transform 



R n — Rn, ®Rnoi Rm = 



5, A j>-l)(f-l/2)K 
— -cos ' 



(2.2) 

J s,t=l 



with 8 s . t denoting the Kronecker symbol. 

The explicit matrix structure is A n = Toeplitz(V) + Hankel((r(V),./(7(V)), with V = 
[VqVi ... V qi 0... 0] and where each V^, i\ = 1 , . . . , qi is the unilevel DCT-III matrix 
associated to the if row of the PSF mask, i.e., Vjj = Toeplitzfv,, ) +Hankel(a(vi 1 ), 
Ja(vi 1 )), with v,-, = [/ifj.o, ■ ■ ■ ,hi u q 2 ,0, . . . ,0]. Here, cr denotes the shift operator such 
thatCT(v fl ) = [Ai 1; 1 , . . . , hi u q 2 , 0, . . . , 0] and J denotes the usual flip matrix; at the block 
level the same operations are intended in block-wise sense. 

Not only the structural characterization, but also the spectral description is com- 
pletely known: let / be the bivariate generating function associated to the PSF mask 
( fT3] l, that is 

9\ 92 

f{x u x 2 ) = hofl+2 £ h Su ocos(sixi) +2 £ h . S2 cos (52^2) 

Jl=l «2=1 
91 92 

+ 4 E L h si,s 2 cos(sixi)cos(s2X2), (2.3) 

S[ = 1 So = l 

then the eigenvalues of the corresponding matrix A„ G %i are given by 

,„] _ (r-\)n 



X s (An)=f(x [ P 1 1 \xf 2 2] ),S=(s l) S 2 ), 
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where s\ — 1,... ,n\, s 2 = 1, . . . ,n 2 , and where the two-index notation highlights the 
tensorial structure of the corresponding eigenvectors. Finally, note that standard op- 
erations like matrix-vector products, resolution of linear systems and eigenvalues 
evaluations can be performed by means of FCT-III flTl within O(nin2log(nin 2 )) 
arithmetic operations (ops). 

2.2 Anti-Reflective boundary conditions 

More recently, Anti-Reflective boundary conditions (AR-BCs) have been proposed 
in 1 20 1 and studied [1.2 3,810.918. 22 1 . The improvement relies in the fact that not 
only the continuity of the image, but also of the normal derivative, are guaranteed 
at the boundary. This higher regularity, not shared with Dirichlet or periodic BCs, 
and only partially shared with reflective BCs, significantly reduces typical ringing 
artifacts in the restored image. 

The key idea is simply to assume that the scene outside the FOV is the anti- 
reflection of the scene inside the FOV. For example, with a boundary at x\ =0 the 
anti-reflective condition imposes /(— X\,x 2 ) — f{x\,x 2 ) = —(f(xi,x 2 ) — f{x*,x 2 )), 
for any x 2 , where x\ is the center of the one-dimensional anti-reflection, i.e., 

/( -xi , x 2 ) = 2/ (x* , x 2 ) - f (jci , x 2 ) , f or any x 2 . 

Notice that, in order to preserve a tensorial structure, a double anti-reflection, first 
with respect to one axis and after with respect to the other, is considered at the corners, 
so that the BCs impose 

f(-Xi , -x 2 ) = 4f(x\,x* 2 ) - 2f(x*,x 2 ) - 2/(*i ,x* 2 ) +f(xi,x 2 ), 

where (x\,x 2 ) is the center of the two-dimensional anti-reflection. More specifically, 
by choosing as center of the anti-reflection the first available data, along the borders, 
the BCs impose 

/l— (l,(2 = 2 1 /l ! ( 2 — /f'i + l,i 2 i /«i+!'i,i2 = 2/« 1 ,i 2 — fn\—i\,ii i*l — l)---i?li h — li---i n 2i 
Ji\ ,1— fe = 2yij,i — fi i ,12+1 1 //i,n2+!2 = ^/'ii«2 — JhAi—h ''1 = li---i w li *2 — lj---i?2- 

At the corners, for any ;'i = 1 , . . . , q\ and i 2 = 1 , . . . , q 2 , we find 

/l-<i,i-/ 2 = 4/i,i -2/ii 2+ i -2/, 1 + i,i +/, 1 +i,, 2 +i, 

f\— il,«2+!2 fl> n 2 ^fi,l2— h ^/il + l,"2 *2' 

fni+i\,l— 12 ^fn\.\ ^fn\ ,(2 + 1 ^ffi[—i[,\ "F fn\ — i\ ,i 2 +l i 
/ni+i'l, 112+12 — ^Jn\,ri2 ^f n l> n l~h ^fn\— ii,n 2 "F /«i — i\ ,« 2 — h " 

As a matter of fact, the rectangular matrix A is reduced to a square Toeplitz-plus- 
Hankel block matrix with Toeplitz-plus-Hankel blocks, plus an additional structured 
low rank matrix. More details on this structure in the general case are reported in 
Section[5] Hereafter, we observe that again under the assumption of strong symmetry 
of the PSF and of a mild finite support condition (more precisely /i, = if \i> \ > n — 2, 
for some j G {1,2}), the linear system A n f = g is such that A„ belongs to the si 
commutative matrix algebra ]2l . 
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This new algebra shares some properties with the T (or DST-I) algebra [5 1. Going 



inside the definition, a matrix A n g 



2D 



has the following block structure 



Ay, 



£o+Zi 


T 





H l +Z 2 







H qi -i +Z qi 


t(H ,...,H qi ) 














H qi -i +Z qi 







H l +Z 2 





6 1 ' 


Ho+Zi 



(2.4) 



where t(Hq, . . . , H qi ) is a block T matrix with respect to the &/ M lD blocks , i\ 
I ..... f/ 1 and = 2Y^L k Ht for k — \ ....,q\. In particular, the s#S% XD block is 



associated to /'/' row of the PSF, i.e., h, 



' 1D ' = [Ai'i,( 2 ]i2=-?2,— m an< ^ ^ is defined as 










hi u l +Zi u 2 







"'1,92-1 + Z <1,92 







,92 




"<1 ,92 


" 




^'1,92-1 + Z l'l.?2 







A<1,1 +Z/i,2 





o 7 ' 


A<1, 0+2/i,l 



(2.5) 



where Zi u k = 2 /i;, ,r for £ = 1, • ■ ■ ,<72 and x(hi u o, . . .,hj im ) is the unilevel T ma- 
trix associated to the one-dimensional PSF /zj| D ' previously defined. Notice that the 
rank-1 correction given by the elements z^k pertains to the contribution of the anti- 
reflection centers with respect to the vertical borders, while the low rank correction 
given by the matrices Z# pertains to the contribution of the anti-reflection centers with 
respect to the horizontal borders. 

Favorable computational properties are guaranteed also by virtue of the z struc- 
ture, so that, firstly we briefly summarize the relevant properties of the two-level T 
algebra [5|. Let ,% = {A„ € R*M xN ^ , n = (m,n 2 ),^(«) = «i«2 | K = Q n A n Q n } 
be the two-level t matrix algebra, i.e., the algebra of matrices that are simultaneously 
diagonalized by the symmetric orthogonal transform 



Qn = Qn, ®Qn v Qm = 



sin 



stn 



(2.6) 



With the same notation as the DCT-III algebra case, the explicit structure of the matrix 
is two level Toeplitz-plus-Hankel. More precisely, A„ =Toeplitz(V) — Hankel((7 2 (V), 
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Jo 2 (V)) with V =[V Vi ... V m 0...0], where each V h , h = is a the 

unilevel x matrix associated to the row of the PSF mask, i.e., Vjj =Toeplitz(v,- 1 ) — 
Hankel((J 2 (v !l ),7(7 2 (vi 1 )) with v i{ = [/!,,, o, ■ ■ . ,/i,- li92 ,0, . . . ,0]. Here, we denote by a 2 
the double shift operator such that (J 2 (v, | ) = [/z^.2, ■ • • ,/i, li92 ,0, . . . ,0]; at the block 
level the same operations are intended in block-wise sense. The spectral character- 
ization is also completely known since for any A„ 6 ST n the related eigenvalues are 
given by 

X s (An)=f(xt[ l] ,X^ ] ),S=(s h S 2 ), xY- 



rK 



l 



where s\ = 1,. .. S2 = I,.. . ,«2, and / is the bivariate generating function associ- 
ated to the PSF defined in ( |2J] >. 

As in the DCT-III case, standard operations like matrix-vector products, resolu- 
tion of linear systems and eigenvalues evaluations can be performed by means of 
FST-I within (9(«i«2log(«i«2)) (ops). Now, with respect to the si ' 0? 2 n D matrix alge- 
bra, a complete spectral characterization is given in [2 3 1. Of considerable importance 
is the existence of a transform T„ that simultaneously diagonalizes all the matrices be- 
longing to £/£% 2D , although the orthogonality property is partially lost. 

Theorem 2.1 Any matrix A n € 2% 2 n D , n = {n\,ri2), can be diagonalized by T n , 
i.e., 

An T n A n T n , T n T n 
where T n =T ni ®T„ 2 , T„ = T ni g3 T„ 2 , with 

r 



a„ 



a m V Qm-2 a m l Jp 



1 



and T m 



~Qm-2P Qm-2 ~Qm-2Jp 







o 7 



The entries of the vector p 6 IR m ~ 2 are defined as p j = 1 — j / (m — 1), j = 1 , . . . , m — 2, 
J £ ]g™-2xm-2 j s t he flip ma trix, and (X m is a normalizing factor chosen such that the 
Euclidean norm of the first and last column ofT m will be equal to 1. 



Theorem 2.2 $2&LetA n e 
[hi 



2D 



n — («i , n2), the matrix related to the PSF hp$p - 



li'2J<l=-9l,-,9li'2=-?2vi92 



. Then, the eigenvalues of A n are given by 



— 1 with algebraic multiplicity A, 

- the «2 — 2 eigenvalues of the unilevel X matrix related to the one -dimensional PSF 
h^''} = Ef 1 = __ hi l ^q 2 , . . . ,Y,V=— a , h-hm]' ea ch one with algebraic multiplicity 2, 



-91 

ll=-<7i '-'l,-<?2' ■ ■ ■ jZ-,' |= - 9i "li.</2J 

- the n i — 2 eigenvalues of the unilevel X matrix related to the one -dimensional PSF 



h {c} = Ef 2 „ h 



-92 
-12 = 



h ■ 1 

-q 2 



each one with algebraic multiplicity 2, 
— the (n\ — 2~){ni — 2) eigenvalues of the two-level X matrix related to the two- 
dimensional PSF hpsp. 



It's worthwhile noticing that the three sets of multiple eigenvalues are related 
to the type of low rank correction imposed by the BCs through the centers of the 
anti -reflections. More precisely, the eigenvalues of X, n -2{h^) and of x, n ^2{h^) 
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take into account the condensed PSF information considered along the horizontal 
and vertical borders respectively, while the eigenvalue equal to 1 takes into account 
the condensed information of the whole PSF at the four corners. In addition, the 
spectral characterization can be completely described again in terms of the generating 



function associated to the PSF defined in (2.3 1, simply by extending to the standard 
T evaluation grid, i.e., it holds 



A S (A„) =f(x^ 1 \x^ y j ,s= {si,s 2 ),Sj = 0,...,rij, 



1 



where the 0— index refers to the first/last columns of the matrix T m 0. See (HO 
for some algorithms related to standard operations like matrix-vector products, res- 
olution of linear systems and eigenvalues evaluations with a computational cost of 
0(n\n2\og{n\ri2)) ops. In fact, the computational cost of the inverse transform is 
comparable with the direct transform one and the very true penalty seems to be the 
loss of orthogonality due to the first/last column of the matrix T m . 

We stress that the latter complete spectral characterization and the related fast 
algorithms for computing the eigenvalues are essential for the fast implementation of 
the regularization algorithms used in the numerical section. 



3 Optimal preconditioning 

In this section we consider in more detail the matrices arising when Anti-Reflective 
BCs are applied in the case of a non-symmetric PSF, the aim being to define the cor- 
responding optimal preconditioner in the srf StiP® algebra. 

More precisely, let A = A(h) be the Anti-Reflective matrix generated by the generic 
PSF h PSF = [hh,h\ h =- qi ,..., quh =- q2 ,... m and let p = p (s) G s^Mf be the Anti- 
Reflective matrix generated by the symmetrized PSF spsf = [ s i l .i 2 ]i l =-q l qi i 2 =-q 2 q 2 
We are looking for the optimal preconditioner P* = P* (s* ) in the sense that 

P* = arg min||A-P||^ , * = argmin||A(/) -P(s)||^ , (3.1) 
where is the Frobenius norm, defined as \\A\\j? = Y. | a !,;'| 2 - Indeed, an anal- 

V'j 

ogous result is know in [17| with respect to Reflective BCs: given a generic PSF 
hpsF — \hi\,h\< tne optimal preconditioner in the DCT-III matrix algebra is generated 
by the strongly symmetric PSF spsf = [<Sii,i 2 ]> given by 

s ±i\.±i 2 = " ^ _ - 

Our interest is clearly motivated by the computational facilities proper of stfSft 1 ® 
algebra, coupled with its better approximation properties. We preliminary consider 
the one-dimensional case in order to introduce the key idea in the proof with a simpler 
notation. Moreover, the proof argument of the two-dimensional case is also strongly 
connected to the one-dimensional one. 
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3.1 One-dimensional case 



2.2 



the 



Let us consider a generic PSF hpsp = [^j],- = _g q - As introduced in Section ! 
idea is to apply an anti-reflection with respect to the border points f\ and /„. Thus, 
we impose 

fi-i = 2/i - /l+i, fn+i = 2/„ - /„-,-, i = 1, . . . , q. 



The resulting matrix shows a more involved structure with respect to the Reflective 
BCs, i.e., it is Toeplitz + Hankel plus a structured low rank correction matrix, as 
follows 



A = 



vo 


T 

u 


' 


Vl 






Vq 


B 


W q 






W\ 





-{Juf 


W 



(3.3) 



with 



u T = [h-i-hi,...,h- q -h q ,Q,...,0], - (Juf = [0,...,0,h q -h^ q ,...,hi -/t_i], 

v k =h k + 2 hj, w k =h_ k + 2 Y, h -j, 
j=k+i j=k+i 

B = T([h- q , . . .,h q ]) -H TL ([h 2 , . . . ,h q \) -H B R([h- 2 , ■ ■ -,h- q ]), 

where T([h- q ,. . . , h q ] ) is the Toeplitz matrix associated to the PSF hpsF, while ( [hi , 
. . . ,h q ]) and //br([/?-2, ■ ■ ■ ,h- q }) are respectively the top-left Hankel and the bottom- 
right Hankel matrices 



hi h3 
h 3 



h q 



hg 



h q 



h q 



••• 



HbR — 



h. 



h- c 



h- 



••• Oh. 



h-3 
h-i h-2 



On the other hand, the Anti-Reflective matrix generated by a strongly 

symmetric PSF spsf = [s q , ■ ■ ■ ,s\,sq,si, . . . ,s q ], among which the minimizer P* in 
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([XT} will be searched, is clearly given by 







Q T 







n 






p = 




t(s) 










n 







6 1 ' 


ro 



where r^ = + 2 £ Si and t(s) is the T (or DST-I) matrix generated by the PSF 

j=k+\ 

SPSF- 

The optimality of the Anti-Reflective matrix generated by the symmetrized PSF 
defined as 

h -i+ h i nA s 

s±i= — ^ ' ^ 3 ' 

can be proved analogously as in ifTTll with respect to the internal part Cj~B t(s) and 
by invoking a non-overlapping splitting argument in order to deal with the external 
border Cb- In fact, we have 



IICII 



\\Cb 



2 

I J? 



and it easy to show that the minimizer found for the first term is the same than for the 
second one. 

Notice that 1 1 u T \ \ 2 and 1 1 — (Ju) T \ \ 2 are constant terms in the minimization process. 
So, as naturally expected, the obtained minimum value will be greater, the greater is 



the loss of symmetry in the PSF. Moreover, with the choice (3.4i, the first and last 
column in Cb share the same norm, i.e., again the most favourable situation. It is worth 
stressing that the minimization process of the second term Cb allows to highlight as 
the tuning of each minimization parameter can be performed just by considering two 
proper corresponding entries in the matrix, i.e., 



(r p -v p ) + (r p -w p ) =0, p = 0, 



where v p and w p are linear combination of the same coefficients with positive and 
negative indices, respectively. Taking this fact in mind, we can now consider a more 
geometrical approach to the proof, that allows to greatly simplify also the proof with 
respect to the minimization of the internal part and can be applied to any type of BCs 
based on the fact that the values of unknowns outside the FOV are fixed or are defined 
as linear combinations of the unknowns inside the FOV. 



Theorem 3.1 Let A 

PSF h PSF = [h,l 



A(h) be the Anti-Reflective matrix generated by the generic 
. The optimal preconditioner in the s^M]^ algebra is the 



matrix associated with the symmetrized PSF sp$F 

h-i+hi 



,Sl,So,Sl,...,Sq], with 



(3.5) 
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Fig. 3.1 A point R, its swapped point R s , the optimal approximation of both Q* . 



Proof Preliminarily, as shown in Figure 3.1 we simply observe that if we consider 
in the Cartesian plane a point R = (R x ,R y ), its optimal approximation Q*, among the 
points Q — (Q x ,Qy) such that Q x = Q y , is obtained as the intersection between the 
line y = x, with the perpendicular line that pass through R, that is 

y-R y = -(x-R x ) 

hence Q* = Q* = (R x +R y ) /2. The same holds true if we consider the swapped point 
R s = (R y ,R x ), since they share the same distance, i.e., d(R,Q*) = d(R s ,Q*). 

Clearly, due to linearity of obtained expression, this result can be extended also 
in the case of any linear combination of coordinates. Thus, by explicitly exploiting 
the structure of A and P, we define as x-coordinate of a point the entry with negative 
index and as y-coordinate of the same point the corresponding entry with positive 
index. For the sake of simplicity we report an example for q = 3, in which we put in 
evidence the x or y coordinate definition, 



C =A-P= 



co{ Q q 0| 0| 

a>l £ 0o 6( q 0f 

ag 0f 9f Bo 9f 0| <of 

' i)\ o: el do c? o)i 

3 °2 Q. Co ®1 



V 3 V 2 V l ®0 



% 

t% Ci v 0o 0f % 0( 

0% 0| §1 O 0f 0f 8% 

e 3 y ej ej O g «f 

03 2 Co ®f 

d)n 



Here, we set the points 

e i = (df,df) = (h- i ,h i ) 

Oi = (eaf,a^) = (h- i + 2 f A_ ; ,/i i +2 £ hj) = (d?+2 £ 0>.o; • 2 £ 0j) 



./='■+! 



;=i+l 



;=i+ 1 



;=i+l 



M = (vf , v?) - (/»_, - a,, k - h- t ) = (ef-e? x ,ef-e° y ) 



and 



Zo 
Zi 

z 2 



(Co*, Co) 

(M) 



(*o- 
(A-i 



h- 2 ,ho- 
— /i 3 , /ii 

-/i 3 ,/zi - 



*2) = 
-A 3 ) 

A-3) 



(00 _ 02 i 0() _ 02 
= (0f _ 03'^1 " 
= (ef-flf ,0? 



03 ) 
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related to the Hankel corrections. The points 0,, Qi, Zi related to the matrix P are de- 
fined in a similar manner, taking into account the strong symmetry property, i.e. they 
have the same x and y coordinates. More in general, the key idea is to transform the 
original minimization problem in the equivalent problem of minimizing the quantity 



c o d(0 o ,0o) 2 + • • • + c q d(0 q ,0 q ) 2 + d(Z ,Z ) 2 + . . . + d(Z m ,Z m ) 2 

+d{n ,n ) 2 + . . . +d(n q ,n q ) 2 +a(N u o) 2 + . . .+d{N q ,o) 2 1 



where Cj are some constants taking into account the number of constant Toeplitz 
entries. Now, by referring to the initial geometrical observation, we start from points 
pertaining to the Toeplitz part, that can be minimized separately, and we obtain the 



minimizer (3.5 I. It is also an easy check to prove the same claim with respect to any 



other terms, by invoking the quoted linearity argument. 



3.2 Two-dimensional case 



be a generic PSF. As introduced in Section 



Let hpsF = \hi, ;J • 

2.2 the idea is to apply an anti -reflection with respect to the border points /i i 2 , fi u \ 
and f ni ,h, fi lt n 2 < h = !>••• 7 n i> h = ,112, and a double anti-reflection at the cor- 
ners in order to preserve the tensorial structure. The resulting matrix shows a more 
involved structure, i.e., it is block Toeplitz + Hankel with Toeplitz + Hankel blocks 
plus a structured low rank correction matrix, as follows 





u 


' 


Vi 






v 91 


B 


w gi 











-JU 


W 



(3.7) 



with 



u= -^,o,...,o], -JU=[o,...,o,ft gi -fi- qi ,- -A-fi-i 



'I] 



Vj =Hj + 2 £ S it W } = H-j + 2 £ J?_,, 

i=j+i '=7+1 

B = T(H- qi ,...,H qi ) -H TL (H2, ■■■,H qi ) -// B r(^-2, ■ • .,H- qi ). 



where T indicates the block Toeplitz matrix, while Htl and Hbr are respectively the 
top-left block Hankel matrix and the bottom-right block Hankel matrix as just previ- 
ously depicted in the unilevel setting and where the block Hj is defined, according to 
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d33), as 





v j,o 


T 
U J 


" 































-(JujY 


Wjfi . 



(3.8) 



with Bj = T(hj- 



92' ■ 



■ > hj,q 2 ) 



Hjl (h j,2,---,h j m ) - H B r (hj-2,--.,hj- 



'12 > 



Refer to (|2.4|i and (|2.5|l for the structure of the matrix P related to a strongly 



symmetric PSF in which the minimizerP*, see (3.1 1, will be searched 



Theorem 3.2 Let A = A{h) be the Anti-Reflective matrix generated by the generic 

PSFhpsp = [Kh]i l =-q u ...,q l ,i 2 =-q 2 ,...,q 2 - 

The optimal pre conditioner in the ss? Si 1 ® algebra is the matrix associated with the 
symmetrized PSF spsf = [si [ i 2 ] i = n „,,-„=_„„ „,> with 



J »l=— «1 *l.»2=— 92>..-.«2 

/i ,—(9 ~i~ ii ~t~ ,— 



l,-'2 



(3.9) 



Proof The proof can be done by extending the geometrical approach just consid- 
ered in the one-dimensional case: we simply observe that if we consider in the 4- 
dimensional space a points = (R x ,R y ,R z ,R w ), its optimal approximation Q* among 
the points Q — (Q x ,Qy,Q z ,Qw) belonging to the line Jz? 



is obtained by minimizing the distance 
d 2 (i?,fl) = 4t 2 -2t{R x +R y 



-R 7 + R» 



Rl 



-Rt+Rt 



This is a trinomial of the form at 2 + j3? + 7, with a > and we find the minimum by 
using the formula for computing the abscissa of the vertex of a parabola 



R x +R y +R z +R H 



f = _JL_ 

2a 4 

Hence the point Q* is defined as Q* = Q* = Q* = Q* w = t* . The same holds true if 
we consider any swapped point R s , not unique but depending on the permutation at 
hand, since they share the same distance, i.e., d(R,Q*) = d(R s ,Q*). Again, thanks 
to the linearity of obtained expression, this result can be extended also in the case of 
any linear combination of coordinates. 
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Thus, by explicitly exploiting the structure of the matrices A and P, we define 
a point by referring to the entry with positive and negative two-index. For instance, 
points pertaining to the Toeplitz part are defined as 

®>i>»2 = ( ^'1,(2' *i ,*2 ' h,h' hfo) ~ { s —h—hi s ~'uh> s ht~'2' s hth)> 
respectively. 

As in the unilevel setting, the original minimization problem is transformed in 
the equivalent problem of minimizing the sum of squared distances analogously as in 
(|3.6|l. We start again from points pertaining to the Toeplitz part, that can be minimized 



separately, and we obtain the minimizer (3.9 1. It is also an easy check to prove the 



same claim with respect to any other couple of points pertaining to Hankel or low 
rank corrections, by invoking the quoted linearity argument. 

It is worth stressing that this proof idea is very powerful in its generality. It can 
be applied to any type of BCs based on the fact that the values of unknowns outside 
the FOV are defined as linear combinations of the unknowns inside the FOV, so that 
it may be useful in the future to prove theoretical results for new proposed BCs. 



4 Computational results 

A well-known iterative method for solving the image deblurring problem is Landwe- 
ber method [ 15 1, whose (k + 1 )-th iteration step is defined by 

x k +i =x k + TA H (g-Ax k ) 1 (4.1) 

where A is the blurring matrix, g is the recorded image and T is the descent parameter 
(we set it equal to one). As one can observe experimentally, the restorations seem to 
converge in the initial iterations, before they become worse and finally diverge; this 
phenomenon is called semiconvergence. Hence Landweber method is a regularization 
method, where the number of steps k is the regularization parameter. Moreover it has 
good stability properties, but it is usually very slow to converge to the sought solution. 
Therefore it is a good candidate for testing the proposed preconditioning technique. 
Thus we introduce the preconditioned Landweber method 

x k + i=x k + XDA H {g -Ax k ), (4.2) 

where D is the preconditioned In order to build it, we compute the eigenvalues Xj 
of the blurring matrix associated to the PSF and to periodic BCs (via FFT) or to the 
symmetrized PSF and to Reflective BCs (via FCT) or to the symmetrized PSF and 
to Anti-Reflective BCs (via FST, see Theorem |2.1| and Theorem |2.2| and comments 
below), then we apply the Tikhonov Filter 
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to determine the eigenvalues dj of D; finally the PSF related to D can be obtained via 
IFFT or IFCT or IFST (the inverses of the previous transforms, namely inverse FFT, 
inverse FCT, inverse FST)). In numerical experiments we have set the parameter a 
manually, so that we have reached excellent performances both in terms of quality of 
the restorations and in acceleration of the method. 

Actually in our implementation, which is partially based on the Matlab Toolbox 
RestoreTools lfT6ll . we have never worked with A H , but always with A', that is the 
matrix related to the PSF rotated by 180 degrees. This approach is known in litera- 
ture as reblurring strategy [10|. The reason behind this choice resides in one of the 
main problems of Anti-Reflective algebra srfffl, i.e. the fact that it is not closed under 
transposition. We stress that A H and A' are the same thing in case of periodic and zero 
boundary conditions, but they are different for Reflective and Anti-Reflective ones. 

To test these different BCs and preconditioning techniques, we have taken into 



account the Cameraman deblurring problem of Figure 4.1 in which the PSF is a 



slightly non-symmetric portion of a Gaussian blur, and the Bridge deblurring problem 
of Figure 4.4 in which the PSF is an highly non-symmetric portion of a Gaussian blur. 
In both cases we have generated the blurred and noisy data g, adding about 0.1% 
of white Gaussian noise. We have chosen to add a low level of noise to emphasize 
the importance of boundary conditions, which play a leading role when the noise 
is low, while they become less decisive when it grows up. Since we know the true 
image /, to measure the quality of the deblurred images we compute the Relative 
Restoration Error (RRE) |]jc — / |j/|| ^, where is the Frobenius norm andx 
is the computed restoration. 

As we expected, from Table |4.l| and Table |4.2] we can notice that both Reflective 
and Anti-Reflective boundary conditions outperform periodic ones, which give rise to 
poor restorations (see first image of Figure 4.2 and Figure 4.5 i. Furthermore by means 
of Anti-Reflective BCs (see third image of Figure 4.2 and Figure 4.5 i we can gain 
restorations of better quality compared with ones obtained employing Reflective BCs 
(see second image of Figure |4~2| and Figure |4~5] >. From Tables [4~T]|4.2| and Figures 
14.314.61 we can see that all these considerations hold also for D-Landweber method 
— i.e. Landweber method with preconditioning — which for a suitable choice of the 
parameter a is able to reach restorations of the same quality of the classical Landwe- 
ber method in much smaller number of steps. In particular the reduction in steps 
for both Reflective and Anti-Reflective BCs is around 50 times for the Cameraman 
deblurring problem and around 8 for the Bridge deblurring problem. 

We stress that the iteration count reported in Table 4.2 in the Anti-Reflective row 
does not have to deceive, because, as it can be seen from Figure 4.7 if we compare 
the restorations gained by Landweber (preconditioned or not) at any given fixed iter- 
ation, employing Reflective BCs or Anti-Reflective BCs, we see that the latter shows 
always equal or better restoration quality. The same remark holds for the Cameraman 



deblurring problem (see Table 4.1 1. In fact Figure 4.7 is very instructive because it 
tells to the generic user two things: a) the curves for Reflective and Anti-Reflective 
BCs are very flat, b) the approximation obtained when using Anti-Reflective BCs is 
always better or equal to that obtained with Reflective BCs. The combined message 
of the previous two items is that we can safely choose the Anti-Reflective BCs, even 
when we are unable to estimate precisely the stopping criterion for deciding the op- 
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timal iteration: we notice that this observation does not hold for the periodic BCs 
where a small error in the evaluation of the optimal iteration leads to a substantial 
deterioration of the quality of the resulting restored image. 

In the end, from the results reported in this section we can say that our proposal 
of the optimal preconditioner in the context of Anti-Reflective BCs is as effective as 
the one introduced in [ 17] for Reflective BCs. Therefore the present work represents 
a theoretical and numerical continuation and strengthening of that line of research. 
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0.1611 
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Anti-Reflective 
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0.1582 
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Table 4.1 Results of the classical and preconditioned Landweber method related to the Cameraman de- 
blurring problem, employing different BCs. 
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Fig. 4.3 Preconditioned Landweber restorations, employing periodic, reflective, anti-reflective BCs. 




Fig. 4.4 Bridge deblurring problem: true image, PSF, blurred and noisy image. 
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Table 4.2 Results of the classical and preconditioned Landweber method related to the Bridge deblurring 
problem, employing different BCs. 




Fig. 4.5 Landweber restorations, employing periodic, reflective, anti-reflective BCs. 

5 Conclusions and Perspectives 

Inspired by the theoretical results on optimal preconditioning stated in [17] in the 
Reflective BCs environment, in this paper we have presented analogous results for 
Anti-Reflective BCs. In both cases the optimal preconditioner is the blurring matrix 
associated to the symmetrized PSF. We stress that our proof is based on a geometrical 




Pietro Dell'Acqua et al. 




Fig. 4.6 Preconditioned Landweber restorations, employing periodic, reflective, anti-reflective BCs. 
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Fig. 4.7 Bridge deblurring problem: RRE trends of Landweber (on the left) and Z)-Landweber (on the 
right) for different BCs. 



idea, which allows to greatly simplify the used arguments, even when non orthogonal 
transforms are involved. Moreover that idea is very powerful in its generality and it 
may be useful in the future to prove theoretical results for new BCs. 

Computational results have shown that the proposed preconditioning strategy is 
effective and it is able to give rise to a meaningful acceleration both for slightly and 
highly non-symmetric PSFs. On the other hand, symmetrization is efficient when we 
have a PSF that is near to be symmetric and it becomes more and more ineffective as 
the PSF departs from symmetry. In this case, other techniques [7 1, which can manage 
directly non-symmetric structures, can gain better performances and in this direction 
we see a substantial development in the near future. 
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